On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments

From: Carrassi, A. et al. Data assimilation in the geosciences: An overview of methods, issues, and perspectives. Wiley Interdisciplinary Reviews: Climate Change 9.5 (2018): e535.

• Numerical weather prediction (NWP) uses filtering techniques to produce weather forecasts in an observation-analysis-forecast cycle.
• In the left, we assume that we have an initial forecast-prior for the physics-based numerical model state.
• This can represent, e.g., the forecasted state of the atmosphere in terms of temperature, pressure and humidity over Reno, six hours in the future.
• This is generated from the simulation of the geophysical fluid dynamics, as they evolve in time given the appropriate forcings and boundary conditions.
• Mathematically, this is described by a system of PDEs, discretized spatially over three dimensions of the atmosphere and surface topography.

From: Carrassi, A. et al. Data assimilation in the geosciences: An overview of methods, issues, and perspectives. Wiley Interdisciplinary Reviews: Climate Change 9.5 (2018): e535.

• We suppose that the “true” physical process evolves in time to the end of the first forecast.
• This corresponds to our real atmosphere reaching the forecast state of plus six hours.
• We suppose that at this time, we receive an observation of the atmosphere’s state variables.

From: Carrassi, A. et al. Data assimilation in the geosciences: An overview of methods, issues, and perspectives. Wiley Interdisciplinary Reviews: Climate Change 9.5 (2018): e535.

• The observation comes with a likelihood function;
  • with the likelihood function we compute the likelihood of the model forecast given the observation.
• Using the forecast-prior and the likelihood, we estimate the Bayesian update of the prior to the posterior conditioned on the observation.
• The posterior density is denoted the analysis.

From: Carrassi, A. et al. Data assimilation in the geosciences: An overview of methods, issues, and perspectives. Wiley Interdisciplinary Reviews: Climate Change 9.5 (2018): e535.

• The analysis is then used as the initialization of the next forecast.
• Forecasting the model state forward in time, we produce a new prior.
• Chaotic evolution (the butterfly effect) that is typical in geophysical problems causes the model forecast to drift from the true state.
• This is in addition to the errors in:
1. accurately representing the true physical process with the imperfect numerical model; and
2. estimating the “true” initial condition from incomplete and noisy observations.

From: Carrassi, A. et al. Data assimilation in the geosciences: An overview of methods, issues, and perspectives. Wiley Interdisciplinary Reviews: Climate Change 9.5 (2018): e535.

• When a new observation becomes available…

From: Carrassi, A. et al. Data assimilation in the geosciences: An overview of methods, issues, and perspectives. Wiley Interdisciplinary Reviews: Climate Change 9.5 (2018): e535.

• We produce a new analysis…

From: Carrassi, A. et al. Data assimilation in the geosciences: An overview of methods, issues, and perspectives. Wiley Interdisciplinary Reviews: Climate Change 9.5 (2018): e535.

• And the cycle continues ad-infinitum…
• When the time series of observations is artificially generated by a numerical model in parallel with the numerical model forecast, this is known as a twin experiment.
• In a twin experiment, a “truth-twin” is run to:
  1. generate the artificial timeseries of observations; and
  2. then to verify the accuracy of the filter’s predictions and analysis.
• The “model-twin” refers to the simulated forecast used in the observation-analysis-forecast cycle.
• Typically, the accuracy of the observation-analysis-forecast cycle is described in the root-mean-square-error (RMSE) of the estimated mean state of the analysis posterior density.
• This is computed using the difference of the model-twin mean state with the known truth-twin state.

Courtesy of: Kaidor via Wikimedia Commons (CC 3.0)

• Imagine a simplified 1-dimensional model of the atmosphere around a lattitude circle.
• We will suppose that this lattitude circle can be discretized into \( n \) total longitude sectors;
• each sector \( i \) is represented by a single state variable \( x_i \).
• The classical Lorenz-96 model is given as, \( \frac{\mathrm{d}\mathbf{x}}{\mathrm{d} t} \triangleq \mathbf{f}(\mathbf{x}), \) where,
• for each state component \( i\in\{1,\cdots,n\} \), \[ \begin{align}\large{f_i(\mathbf{x}) =-x_{i-2}x_{i-1} + x_{i-1}x_{i+1} - x_i + F}.\end{align} \]
• The state variables \( x_i \) have periodic boundary conditions modulo \( n \), \( x_0=x_n \), \( x_{-1}=x_{n-1} \) and \( x_{n+1}=x_{1} \).
• The term \( F \) in the Lorenz-96 system is the forcing parameter that injects energy to the model.
• The terms in the Lorenz-96 model, while not a true physics-based model, approximate geophysical behavior with:
1. external forcing and internal dissipation with the linear terms; and
2. advection and conservation of energy in the quadratic terms.

From: Wilks, D. Effects of stochastic parametrizations in the Lorenz'96 system. Quarterly Journal of the Royal Meteorological Society 131.606 (2005): 389-407.

• A two layer version of the Lorenz-96 system also exists which simulates coupled, ocean-atmosphere dynamics.
• The inner circle represents oceanic variables evolving at a slow time scale;
• these are coupled to the outter circle of atmospheric variables evolving at a fast time scale.
• The two layer model is of particular interest for examining the effects of model uncertainty,
• e.g., when the fast variables are not evolved accurately in the model-twin.
• This is a common model to study stochastic model reduction.
• E.g., instead of simulating the fast-variable dynamics, one may parameterize their effect on the slow variables with a stochastic process.
• This kind of approximation can be justified mathematically due to the central limit theorem:
• in the asymptotic separation of the time scales for the fast and slow variables,
• the effect of the fast variables on the slow variables will reduce to additive, Gaussian noise which perturbs a mean-field ordinary differential equation².

From: Grudzien et al. On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments. In submission.

• In the left, we plot the estimated strong convergence discretization error in the vertical axis versus the time-discretization size in log-log base 10 scale.
• The point estimates are given by the average estimated discretization error over 500 independent batches of simulations.
• Each batch simulates the evolution of an initial condition with respect to 100 independent Wiener processes.
• For each Wiener process, we simulate a finely discretized path solution \( \mathbf{x}_\mathrm{SP} \) with error on \( \mathcal{O}\left(10^{-7}\right) \), and a coarser approximation by one of the tested methods.
• We estimate the quantity, \[ \begin{align} \mathbb{E}\left[\left\Vert \mathbf{x}(T) - \mathbf{x}_\mathrm{SP}(T)\right\Vert\right] \end{align} \] as the mean difference of the finely discretized solution versus the coarse discretization over all Wiener processes in the batch.

• It is known that these batch-mean estimates are realizations of Gaussian random variables distributed around the true expectation for sufficiently many samples of Wiener processes.
• Therefore, the slope (order of convergence) is estimated with weighted least squares with weights proportional to the inverse standard deviation of the batch realizations.

From: Grudzien et al. On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments. In submission.

• The parameter \( s \), fixed for each panel, is the diffusion level in the L96-s model.
• We note that all the orders of convergence are empirically verified here.
• However, the constant \( C \) in the bound \( C\Delta^\gamma \) has a large impact on the overall strong discretization error.
• Particularly, we see that in the low diffusion regime, the order 1.0 Runge-Kutta scheme has an overall discretization error comparable to the order 2.0 Taylor scheme.
• The performance of the Runge-Kutta scheme varies significantly, however, between low and high diffusion levels.

• While the order of convergence remains the same in all diffusion regimes, the difference in performance is reflected in the constant \( C \).

From: Grudzien et al. On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments. In submission.

• The effect of the constant is even more pronounced in the weak convergence metric.
• We see that the order 2.0 Taylor scheme fails to outperform the order 1.0 Runge-Kutta scheme in this metric for the low diffusion regime.
• Moreover, the overall weak discretization error varies on an order of magnitude between low and high diffusion levels.
• The high precision in the strong sense across all diffusion regimes makes the Taylor scheme a natural choice for generating the truth-twin.
• On the other hand, the Runge-Kutta scheme is appealing to generate the ensemble forecast of the model-twin, due to its accuracy in the weak sense.

• However, because the weak convergence of the Runge-Kutta scheme varies by an order of magnitude between diffusion levels, we use the Taylor scheme as a benchmark in the following experiments.
• We study how the ensemble forecast statistics of the Euler-Maruyama and the Runge-Kutta scheme differ from those generated by the Taylor scheme.

From: Grudzien et al. On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments. In submission.

• In the left plot, we generate 500 initial conditions of the L96-s model and forecast the initial condition with 100 independent Wiener processes to time plus 10 (horizontal axis).
• This ensemble of evolved states is then used to compute the empirical, ensemble-estimated mean and ensemble-estimated spread of the forward distribution.
• This is performed over each of the integration schemes:
1. order 2.0 Taylor
2. Euler-Maruyama
3. 4-stage Runge-Kutta
• In the top panels, we compute the root-mean-square-deviation of the ensemble mean of either the Euler-Maruyama scheme or the Runge-Kutta scheme from the ensemble mean of the Taylor scheme as a base of comparison.
• In the bottom panels, we compute the ratio of the ensemble spread of either the Euler-Maruyama scheme or the Runge-Kutta scheme with the ensemble spread of the Taylor scheme as a basis of comparison.
• Here, the Taylor scheme uses a step size of \( \Delta=10^{-3} \) and the same step size is used for each the Euler-Maruyama and Runge-Kutta schemes.

From: Grudzien et al. On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments. In submission.

• We plot the summary statistics versus forecast time over the 500 initial conditions, where:
• the median is plotted as a solid line;
• the inner 80 percentile is plotted shaded;
• the min/ max values are plotted in dashed.
• We note, the ensemble statistics of the Runge-Kutta scheme have almost no appreciable difference from the Taylor scheme for forecast times of up to almost \( T=3 \).
• On the other hand, there is a rapid divergence of the ensemble mean and a bias in the ensemble spread of the Euler-Maruyama, with the median spread ratio above 1.0 asymptotically.

From: Grudzien et al. On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments. In submission.

• On the left, we evaluate the same plot except where the Rung-Kutta and Euler scheme use a step size of \( \Delta=10^{-2} \) while the benchmark Taylor scheme uses a step size of \( \Delta=10^{-3} \).
• In this case, there is a more rapid onset of the divergence of the ensemble means for both schemes, as well as differences induced in the ensemble spread.
• The Euler-Maruyama scheme in particular goes through extreme divergence of the ensemble means, especially in the low difussion regimes.
• Likewise, the ensemble spread is extremely biased for the Euler-Maruyama scheme, where the ensemble spread ratio has a miniumum value consistently above 1.0.
• On the other hand, the Runge-Kutta scheme does not have such extreme divergence of trajectories, and the spread remains unbiased asymptotically.

From: Grudzien et al. On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments. In submission.

• On the left, we plot the twin-experiment RMSE and ensemble spread for an ensemble Kalman filter, in which both the truth-twin and model-twin are generated with the order 2.0 Taylor scheme, each with step size \( \Delta=10^{-3} \).
• The vertical axis is the diffusion level \( s \) while the horizontal axis is the variance of the error in the observations given to the ensemble Kalman filter.
• We take \( N=10^2 \) samples in the ensemble Kalman filter, while the state dimension is \( n=10 \) to guarantee convergence of the filtering statistics without additional techniques (inflation/ localization).
• The RMSE and spread are computed as the asymptotic average over \( 2.5\times 10^{4} \) analysis cycles.

• We see that the analysis RMSE is comparable to the analysis spread, and lower than the standard deviation of the error in the observations.
• This indicates that this is indeed represents the stable filter statistics.
• In the following, we demonstrate the effect of varying the model-twin integration scheme.
• We will vary the method of ensemble integration between the Runge-Kutta scheme and the Euler-Maruyama scheme, each with a step size between \( \Delta\in\left\{10^{-3},10^{-2}\right\} \).

From: Grudzien et al. On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments. In submission.

• As was indicated with the ensemble statistics, the difference between generating the ensemble with the Runge-Kutta with a step size \( \Delta= 10^{-3} \) scheme and the benchmark Taylor scheme is negligible.
• The difference in the filtering statistics for these two configurations is on the order of \( 10^{-6} \).
• We increase the step size of the Runge-Kutta scheme to \( \Delta=10^{-2} \) and plot the differences between the RMSE and the ratio of the ensemble spread, versus the benchmark Taylor configuration.
• When we increase the step size of the Runge-Kutta scheme to \( \Delta=10^{-2} \), we introduce small errors,
• however, the residuals are unstructured in sign or magnitude.

• Indeed, the mean of the difference in the RMSE is approximately \( 8\times 10^{-5} \) though with a standard deviation of approximately \( 10^{-3} \) (on the order of the weak-discretization error).
• This indicates that the observed differences in the filtering statistics can be attributed to random numerical errors, which are themselves unbiased.

From: Grudzien et al. On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments. In submission.

• In contrast, when we generate the ensemble forecast with the Euler-Maruyama scheme with time step \( \Delta=10^{-3} \), we notice structure in the residuals.
• For low diffusion regimes, the RMSE of the ensemble Kalman filter using the Euler-Maruyama scheme suffers from the bias observed in the ensemble forecast.
• Likewise, the spread in the ensemble Kalman filter suffers from this artificial inflation observed earlier.
• However, it is important to note that in high diffusion regimes, the RMSE and spread of the ensemble are nearly identical for the ensemble Kalman filter when using either the Taylor or Euler-Maruyama scheme.

• Particularly, the level at which the bias impacts the twin experiment statistics depends strongly on the observation uncertainty, and most especially on the model uncertainty.
• We investigate next the effect of increasing the step size of the Euler-Maruyama scheme to \( \Delta=10^{-2} \) when generating the ensemble.

From: Grudzien et al. On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments. In submission.

• On the left, we plot the differences in the RMSE and the ensemble spread when generating the ensemble with the Euler-Maruyama scheme with the coarse step-size.
• The structure of the residuals remains largely the same, but the magnitude changes dramatically.
• In particular, for low diffusion regimes we see that the discretization error is sufficient to introduce filter divergence.
• However, the effect of this bias on the filtering statistics is once again relaxed in the presence of high model uncertainty.
• For high levels of diffusion, \( s\geq 0.5 \), the difference of the RMSE with that when the ensemble is generated with the benchmark Taylor scheme is on the order of \( 10^{-2} \).

• However, we consistently see the bias in the ensemble spread, present once again from the earlier ensemble forecast statistics.

From: Grudzien et al. On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments. In submission.

• As an additional check, we test the benchmark configuration versus using:
• the Runge-Kutta scheme with step size \( \Delta=10^{-2} \) to generate the model-twin; and
• the Taylor scheme with step size \( \Delta = 5 \times 10^{-3} \) to generate the truth-twin.
• The configuration ensures that the expected discretization error is less than \( 10^{-3} \) over all diffusion regimes.
• While it is possible that the innacuracies in the truth-twin have now added some structure in the residuals, the results are largely the same as in the configuration with the more accurate truth-twin.

• When increasing the step size of the truth-twin with the model-twin generated by the Euler-Maruyama scheme the results are also largely the same.
• This has the effect of relaxing the issues in the Euler-Maruyama scheme when the model-twin uses a step size of \( \Delta=10^{-3} \).
• However, there is still filter divergence when the model-twin uses a step size of \( \Delta=10^{-2} \).
• These results thus suggest a computationally efficient and statistically robust configuration for twin experiments using the Taylor scheme for the truth-twin and the Runge-Kutta scheme for the model-twin.
• The step size of the truth-twin can be taken as large as \( \Delta=5\times 10^{-3} \) and the step size of the model-twin as large as \( \Delta=10^{-2} \) without biasing the benchmark statistics.